
setwd("~/Desktop/")    
getwd()
##import data
A<- read.delim("A.csv",header=T,sep=',',fill=T)
B<- read.delim("B.csv",header=T,sep=',',fill=T)
C<- read.delim("C.csv",header=T,sep=',',fill=T)
D<- read.delim("D.csv",header=T,sep=',',fill=T)
E<- read.delim("E.csv",header=T,sep=',',fill=T)
F<- read.delim("F.csv",header=T,sep=',',fill=T)


C1<- read.delim("C1.csv",header=T,sep=',',fill=T)
D1<- read.delim("D1.csv",header=T,sep=',',fill=T)
E1<- read.delim("E1.csv",header=T,sep=',',fill=T)
F1<- read.delim("F1.csv",header=T,sep=',',fill=T)

##pooling effect sizes (Random-effects model)
library(meta)
library(metafor)
##A relation
m.A <- metagen(TE,
               seTE,
               data = A,
               studlab = paste(Author),
               comb.fixed = FALSE,
               comb.random = TRUE,
               method.tau = "SJ",
               hakn = TRUE,
               prediction = TRUE,
               sm = "SMD")
m.A
##B relation
m.B <- metagen(TE,
               seTE,
               data = B,
               studlab = paste(Author),
               comb.fixed = FALSE,
               comb.random = TRUE,
               method.tau = "SJ",
               hakn = TRUE,
               prediction = TRUE,
               sm = "SMD")
m.B
##C relation
m.C <- metagen(TE,
               seTE,
               data = C,
               studlab = paste(Author),
               comb.fixed = FALSE,
               comb.random = TRUE,
               method.tau = "SJ",
               hakn = TRUE,
               prediction = TRUE,
               sm = "SMD")
m.C

##D relation
m.D <- metagen(TE,
               seTE,
               data = D,
               studlab = paste(Author),
               comb.fixed = FALSE,
               comb.random = TRUE,
               method.tau = "SJ",
               hakn = TRUE,
               prediction = TRUE,
               sm = "SMD")
m.D

##E relation
m.E <- metagen(TE,
               seTE,
               data = E,
               studlab = paste(Author),
               comb.fixed = FALSE,
               comb.random = TRUE,
               method.tau = "SJ",
               hakn = TRUE,
               prediction = TRUE,
               sm = "SMD")
m.E

##F relation
m.F <- metagen(TE,
               seTE,
               data = F,
               studlab = paste(Author),
               comb.fixed = FALSE,
               comb.random = TRUE,
               method.tau = "SJ",
               hakn = TRUE,
               prediction = TRUE,
               sm = "SMD")
m.F

##meta-regression for A uncoditional model
library(metafor)
model1<- rma.mv(  yi=TE, V=seTE^2,  random = list(~ 1 |ID),
                  data = A)
summary(model1)

##meta-regression for A (whole model)

library(metafor)
modelA<- rma.mv(mods = ~ Study+Country+Disconfirmation+Expectation,
                yi=TE, V=seTE^2,  random = list(~ 1 |ID),
                data = A)
summary(modelA)

(sum(model1$sigma2) - sum(modelA$sigma2)) / sum(model1$sigma2)

library(clubSandwich)
robust.rma.mv(modelA, cluster=A$ID)

##meta-regression for B unconditional

model2<- rma.mv(  yi=TE, V=seTE^2,  random = list(~ 1 |ID),
                  data = B)
summary(model2)
##meta-regression for B (whole model)

modelB<- rma.mv(mods = ~ Study+Country+Disconfirmation+Expectation,
                yi=TE, V=seTE^2,  random = list(~ 1 |ID),
                data = B)
summary(modelB)
(sum(model2$sigma2) - sum(modelB$sigma2)) / sum(model2$sigma2)

library(clubSandwich)
robust.rma.mv(modelB, cluster=B$ID)

##meta-regression for C unconditional

model3<- rma.mv(  yi=TE, V=seTE^2,  random = list(~ 1 |ID),
                  data = C1)
summary(model3)

##meta-regression for C (whole model)

modelC<- rma.mv(mods = ~ Policy+Study+Country+Disconfirmation+Expectation,
                yi=TE, V=seTE^2,  random = list(~ 1 |ID),
                data=C1)
summary(modelC)
(sum(model3$sigma2) - sum(modelC$sigma2)) / sum(model3$sigma2)

library(clubSandwich)
robust.rma.mv(modelC, cluster=C1$ID)

##meta-regression for D unconditional

model4<- rma.mv(  yi=TE, V=seTE^2,  random = list(~ 1 |ID),
                  data = D1)
summary(model4)
##meta-regression for D (whole model)

modelD<- rma.mv(mods = ~Country+Disconfirmation+Expectation,
                yi=TE, V=seTE^2,  random = list(~ 1 |ID),
                data=D1)
summary(modelD)
(sum(model4$sigma2) - sum(modelD$sigma2)) / sum(model4$sigma2)

library(clubSandwich)
robust.rma.mv(modelD, cluster=D1$ID)

##meta-regression for E unconditional

model5<- rma.mv(  yi=TE, V=seTE^2,  random = list(~ 1 |ID),
                  data = E1)
summary(model5)
##meta-regression for E (whole model)

modelE<- rma.mv(mods = ~Study+Country+Disconfirmation+Expectation,
                yi=TE, V=seTE^2,  random = list(~ 1 |ID),
                data=E1)
summary(modelE)
(sum(model5$sigma2) - sum(modelE$sigma2)) / sum(model5$sigma2)

library(clubSandwich)
robust.rma.mv(modelE, cluster=E1$ID)

##meta-regression for F unconditional

model6<- rma.mv(  yi=TE, V=seTE^2,  random = list(~ 1 |ID),
                  data = F1)
summary(model6)
##meta-regression for F (whole model)

modelF<- rma.mv(mods = ~ Study+Policy+Country+Disconfirmation+Expectation,
                yi=TE, V=seTE^2,  random = list(~ 1 |ID),
                data=F1)
summary(modelF)
(sum(model6$sigma2) - sum(modelF$sigma2)) / sum(model6$sigma2)

library(clubSandwich)
robust.rma.mv(modelF, cluster=F1$ID)

##testing for publication bias--funnel plot
funnel(m.A,xlab = "A:Expectation-Disconfirmation Relation 
    Egger's Asymmetry Test, p-value=0.05",ylab="",cex.lab=1.3,mgp = c(3, 0.5, 0),xlim = c(-1.5,1.5),
       ylim = c(0.1,0))
##publication bias for B
funnel(m.B,xlab = "B:Performance-Disconfirmation Relation 
       Egger's Asymmetry Test, p-value=0.07",ylab="",cex.lab=1.3,mgp = c(3, 0.5, 0),xlim = c(-1.5,1.5),
       ylim = c(0.1,0))
##publication bias for C
funnel(m.C,xlab ="C:Disconfirmation-Satisfaction Relation 
       Egger's Asymmetry Test, p-value=0.52",ylab="",cex.lab=1.3,mgp = c(3, 0.5, 0),xlim = c(-1.5,1.5),
       ylim = c(0.1,0))
##publication bias for D
funnel(m.D,xlab = "D:Performance-Expectation Relation  
       Egger's Asymmetry Test, p-value=0.25",ylab="",cex.lab=1.3,mgp = c(3, 0.5, 0),xlim = c(-1.5,1.5),
       ylim = c(0.1,0))
##publication bias for E
funnel(m.E,xlab = "E:Performance-Satisfaction Relation 
       Egger's Asymmetry Test, p-value=0.004",ylab="",cex.lab=1.3,mgp = c(3, 0.5, 0),xlim = c(-1.5,1.5),
       ylim = c(0.1,0))
##publication bias for F
funnel(m.F,xlab = "F:Expectation-Satisfaction Relation 
      Egger's Asymmetry Test,p-value=0.97",ylab="",cex.lab=1.3,mgp = c(3, 0.5, 0),xlim = c(-1.5,1.5),
       ylim = c(0.1,0))

##testing for publication bias--egger's test and SQR of sample size
regA1 = lm(Tese ~ srsample, data = A)
summary(regA1)
regA2 = lm(Tese ~ seinverse, data = A)
summary(regA2)

regB1 = lm(Tese ~ srsample, data = B)
summary(regB1)
regB2 = lm(Tese ~ seinverse, data = B)
summary(regB2)

regC1 = lm(Tese ~ srsample, data = C)
summary(regC1)
regC2 = lm(Tese ~ seinverse, data = C)
summary(regC2)

regD1 = lm(Tese ~ srsample, data = D)
summary(regD1)
regD2 = lm(Tese ~ seinverse, data = D)
summary(regD2)

regE1 = lm(Tese ~ srsample, data = E)
summary(regE1)
regE2 = lm(Tese ~ seinverse, data = E)
summary(regE2)

regF1 = lm(Tese ~ srsample, data = F)
summary(regF1)
regF2 = lm(Tese ~ seinverse, data = F)
summary(regF2)